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Abstract A boundary element method (BEM) simulation is used to compare the efficiency of numerical in- 
verse Laplace transform strategies, considering general requirements of Laplace-space numerical approaches. 
The two-dimensional BEM solution is used to solve the Laplace-transformed diffusion equation, producing 
a time-domain solution after a numerical Laplace transform inversion. Motivated by the needs of numer- 
ical methods posed in Laplace-transformed space, we compare five inverse Laplace transform algorithms 
and discuss implementation techniques to minimize the number of Laplace-space function evaluations. We 
investigate the ability to calculate a sequence of time domain values using the fewest Laplace-space model 
evaluations. We find Fourier-series based inversion algorithms work for common time behaviors, are the 
most robust with respect to free parameters, and allow for straightforward image function evaluation re-use 
across at least a log cycle of time. 

Keywords numerical Laplace transform inversion ■ boundary element method ■ 2D diffusion ■ Helmholtz 
equation ■ Laplace-space numerical methods ■ groundwater modeling 

1 Introduction 

Simulation methods that are posed in Laplace-transformed space, then numerically inverted back to the 

time domain (i.e., Laplace-space methods), are a viable alternative to the more standard use of finite 
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4 differences in time. We use the the two-dimensional boundary element method (BEM) as an example of 

5 this type of approach, to solve the Laplace-transformed diffusion equation (i.e., the Yukawa or modified 

6 Helmholtz equation). We investigate five numerical inverse Laplace transform methods and implementa- 

7 tion approaches, namely the methods of [Stehfest(1970)] , [Schapery(1962)| , [Weeks(1966)] , [Talbot(1979)] , 



a and de Hoog et al(1982) . Naively implemented Laplace-space simulations can be more computationally 

i expensive than using finite differences in time, but they have the advantage of allowing evaluation at 

io any time, without evolving from an initial condition, and image function calculations are trivially par- 

u allelized across Laplace parameters [Davies and Cra nn(2002)|. When Laplace-space numerical models are 

12 used in parameter estimation, hundreds or thousands of forward simulations may be required — making 

o forward model efficiency critical. Although parameter estimation may be done directly in Laplace space 

14 Barnhart and Illangasekare ( 201 2 )] , choosing an efficient inversion strategy is important in most applica- 



tions. 



The Laplace transform has a long history of use to derive analytical solutions to diffusion and wave 
problems (e.g., see list of citations by Duffy(2004) pp. 191-220]). Often the analytical inverse trans- 



form is too difficult to find or evaluate in closed form. A researcher then resorts to approximate ana- 
lytical methods (e.g., [Hantush(1 960),Ste rnberg(1969)] ) or numerical inversion (e.g., [Malama et al(2009)[ 



20 Mishra and Neuman(2010) ). Numerical methods can similarly benefit from the Laplace transform, convert- 

21 ing the time-dependence of a differential equation to parameter dependence. Laplace-space finite-element 

22 approaches have seen application to groundwater flow and solute transport (e.g., |Sudicky and McLaren(1992)| 



Morales-Casique and Neuman(2009) ) , and Laplace-space BEM has also been used in groundwater applica- 



tions (e.g., |Kythe(1995)| §10.3] or |Liggett and Liu(1982) §10.1]). The Laplace transform analytic element 



25 method [Kuhlman and Neum an(2009) is a transient extension of the analytic element method. These dif- 

26 ferent Laplace-space approaches may have diverse spatial solution strategies, but they have a common 

27 requirement of effective Laplace transform numerical inversion algorithms. We couple a BEM model in the 
2a Laplace domain with a numerical Laplace transform inversion routine, but our conclusions should be valid 

29 for both gridded and mesh-free Laplace-space numerical methods. Any Laplace-space numerical approach 

30 begins with determination of optimal Laplace parameter values. Then each image function evaluation is 

31 computed from the simulation. The final step involves approximating the time-domain solution from the 

32 vector of image function values using the algorithm of choice. 



|Bellman et al(1966)] was an early review book on numerical Laplace transform inversion for linear 
and non-linear problems, but without the benefit of the many algorithms that have since been devel- 



oped. Davies and Martin(1979) performed a thorough survey, assessing numerical Laplace transform in- 



version algorithm accuracy for techniques available in 1979, using simple functions for their benchmarks. 



37 Duffy (1993) reviewed the numerical inversion characteristics for more pathological time behaviors using 

3a the Fourier series, Talbot, and Weeks inversion methods. The review book by [Cohen(200f)l summarizes 

39 historical reviews and discusses commonly used inversion their variations. More details and examples can 

40 be found in these reference regarding the convergence behavior of the five inversion algorithms discussed 

41 here. 

42 While these published numerical inverse Laplace transform algorithm reviews are thorough and useful, 

43 they focus on computing a single time-domain solution as accurately as possible. These reviews did not try 

44 to minimize Laplace-space function evaluations, since their functions were simple closed-form expressions, 

45 not simulations. We investigate Laplace transform inversions algorithms that can compute a sequence of 

46 time domain values using the fewest Laplace-space model evaluations possible, a desirable property for use 

47 in Laplace-space numerical methods. Using numerical Laplace transform inversion in a simulation approach, 

48 rather than a time-marching method, allows the researcher to readily switch between fast and accurate by 

49 changing the number of approximation terms in the inversion. 

50 In the next section we define the mathematical formulation of the governing equation and Laplace 

51 transform. In the third section we introduce the five inverse Laplace transform algorithms. In the final 

52 section we compare results using five different inversion algorithms to invert the BEM modified Helmholtz 

53 solution on the same domain with four different boundary conditions, leading to recommendations for 

54 Laplace-space numerical approaches. 



55 2 Governing Equation and Laplace Transform 

56 The BEM model generally simulates substance flow (e.g., energy or groundwater), which can be related to 

57 a potential <f> (e.g., temperature or hydraulic head). The medium property a is diffusivity [L 2 /T], the ratio 

58 of the conductance in the substance flux and potential gradient relation (e.g., Fourier's or Darcy's law) 

59 to the substance capacity per unit mass (e.g., heat capacity or storativity) . The BEM (e.g., [Kythe(1995)[ 
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60 Liggett and Liu(1982),Breb bia et al(1984)| ) can be used to solve the diffusion equation 



61 where a is a real constant in space and time. We consider (QJ in a domain subject to a combination of 

62 Dirichlet <j> (r u (s),t) = f u (s,t) and Neumann n ■ V0(r q (s),t) = / q (s,t) boundary conditions along the 

63 perimeter of the 2D domain r = r u U -T q , where h is the boundary unit normal, and s is a boundary 

64 arc-length parameter. Without loss of generality, we only consider homogeneous initial conditions. 

65 The Laplace transform is 

£{/(*)} = f(p) = / /(t) e - pt dt, (2) 
Jo 

66 where p is the generally complex- valued Laplace parameter, and the over-bar denotes a transformed variable. 

67 The transformed diffusion equation with zero initial conditions is the homogeneous Yukawa or modified 

68 Helmholtz equation, 

V 2 ^-g 2 ^ = 0, (3) 

69 where q 2 = p/a. Equation [3] arises in several groundwater applications, including transient, leaky, and lin- 

70 earized unsaturated flow [Bakker and Kuhlman(2011)| . The transformed boundary conditions are <j> (-Tu(s)) = 

71 /u( s )/t(p) an d n ■ V(f>(rq(s)) = /q(s)/t(f>), where the temporal and spatial behaviors have been decomposed 

72 as in separation of variables. Arbitrary time behavior can be developed through convolution in t (Duhamel's 

73 theorem), which is multiplication of image functions in Laplace space. Here, ft(p) represents the Laplace 

74 transform of the time behavior applied to the boundary conditions. The Laplace transformation makes it 

75 possible to solve transient diffusion (a parabolic equation) using the BEM, which is well-suited for elliptical- 

76 type equations. 

77 The inverse Laplace transform is defined as the Bromwich contour integral, 

CT 1 {/(p) } = f{t) = ± / f(p)e pt dp, (4) 

78 where the abscissa of convergence a > is a real constant chosen to put the contour to the right of all 

79 singularities in f(p). In Laplace-space numerical approaches, ([3} is solved by a suitable numerical method, 
so therefore only samples of f(p) are available; this precludes an analytical inversion. Five numerical inverse 
si Laplace transform algorithms are discussed in the following section. 



82 3 Numerical Inverse Laplace Transform Methods 

S3 Equation [3] is an integral equation for unknown f(t) given /(p); its numerical solution is broadly split into 

84 two categories. Methods are either based on quadrature or functional expansion using analytically invertible 

85 basis functions. [Da vies(2005), Chap. 19] relates most major classes of inverse Laplace transform methods 

86 using a unified theoretical foundation; we adopt a simplified form of their general notation. The Fourier 

87 series and Talbot methods are quadrature-based, directly approximating Weeks' and Piessen's methods 

88 are f(p) expansions using complex-valued basis functions, while the Gaver-Stehfest and Schapery methods 

89 use real-valued functions to accomplish this. 

so The numerical inverse Laplace transform is in general an ill-posed problem (e.g., | Al- Shuaibi (2001)) ). 

91 No single approach is optimal for all circumstances and temporal behaviors, leading to the diversity of 

92 viable numerical approaches in the literature (e.g., [Cohen(2007)| ). 

93 3.1 Gaver-Stehfest Method 



The Post-Widder formula [Widder(194"T)||Al-Shuaibi(2001)] is an approximation to @ that only requires 



95 f(p) f° r rea l P to represent (J2j) as an asymptotic Taylor series expansion. The formula requires high-order 

96 analytic image function derivatives, and is impractical for numerical computation. Stehfest proposed a 

97 discrete version of the Post-Widder formula using finite differences and Salzer summation [Stehfest (1 970)] , 

f(t,N) = ^J2v k f(k^). (5) 
fc=i ^ ' 

98 The Vjt coefficients only depend on the number of expansion terms, TV (which must be even), which are 

min(fe,iV/2) .£/•«, -M 

99 These become very large and alternate in sign for increasing k. The sum ([5]) begins to suffer from cancellation 

100 for TV > the number of decimal digits of precision (e.g., double precision = 16). For ft(p) that are non- 

101 oscillatory and continuous, N < 18 is usually sufficient |Stehfest(1970) . If programmed using arbitrary 



precision (e.g. Mathematica or a multi-precision library [Bailey et al(2002j||Johansson(2011) ), the method 



can be made accurate for most cases [Abate and V alk6(2004)]. Unfortunately, p is explicitly a function of 
t; for each new t, a new f(p) vector is needed. In Laplace-space numerical approaches, each vector element 
is constructed using a simulation, therefore this can be a large penalty. 



106 The method is quite easy to program; the Vj can be computed once and saved as constants. This method 

107 has been popular due to its simplicity and adequacy for exponentially decaying ft (p) . 



3.2 Schapery's Method 



We can expand the deviation of f(t) from steady-state / s using exponential basis functions |Schapery(1962)| , 

N 

f(t,N) = f a +J2^ K \ (7) 
i=i 

where a, is a vector of unknown constants. Applying (|2j) to (JTJ) gives 

(8) 



/(p.,JV) = A+V-5— j = l,2,...,M. 



N 



in The pj are selected (a geometric series is recommended [Ligge tt" and Liu(1982)| ) to cover the important 

112 fluctuations in f(p). After setting pi = pj the aj coefficients can be determined as the solution to PjjOj = 

113 (f(pj) — fs/Pj)- The symmetric matrix to decompose is 

(2 Pl y 1 (pi+psr 1 ... (pi+Pn)- 1 

(P2+Piy 1 (2p 2 ) _1 ... (P2+PN)' 1 
(PN+piy 1 (PN+P2)- 1 ... (2PAT)- 1 

114 which only depends on pj; it can be decomposed independently of f(p) and / s . 

115 This method is not difficult to implement when existing matrix decomposition libraries are available, and 
lie only requires real computation. The method has been used for inverting BEM results |Liggett and Liu(1982)] , 
117 but has two main drawbacks. First, in its formulation above, it requires a steady-state solution, but (0 
us could be posed without / s . Secondly, no theory is presented for optimally picking pj; some trial and error 
119 is required [Liggett and Liu(1982)] . 



120 3.3 Mobius Transformation Methods 

121 We can use the Mobius transformation to conformally map the half-plane right of a to the unit disc, mak- 

122 ing the Laplace domain more amenable to approximation using orthonormal polynomials (e.g., Chebyshev 

123 [Piessens(1972)] , [Lanczos(1988)l §28] or Laguerre [Weeks(1966)] , [Lyness and Giunta(1986)] , [Lanczos(1988)l 

124 §30]). If a was chosen properly, f(p) is guaranteed to be analytic inside the unit circle. The most-used inverse 
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125 Laplace transform method from this class is Weeks' method, which uses a complex power series to expand 

126 f(p) inside the unit circle. Upon inverse Laplace transformation, the power series becomes a Laguerre 

127 polynomial series. 

128 Weeks method is 

N 

/(t,iV + l) = e (K - b / 2)t ^a„L„(6t), (9) 



where L n (z) is an n-order Laguerre polynomial and k and b are free parameters. Weeks suggested k = 
a + 1/tmax and b = ff/t max , where i m ax is the maximum transformed time. The parameters b and k are 
chosen to optimize convergence; some schemes are given [Weidem an(1999) for finding optimum parameter 



values for a given ft{p), but search techniques require hundreds of f(p) evaluations. A more general form of 
([9|) can also be used, which allows for more general asymptotic behavior of the image function [Davies(2005) 



134 §19.5]. Weeks assumed pf(p) is analytic at infinity. The Laplace transform of ([9} is known, but to make it 

135 easier to represent with polynomials, f(p) is mapped inside the unit circle via z = (p — k — 2b) /(p — k + 26) . 

136 The coefficients a n are determined by the midpoint rule, 

1 M_1 

an= 2M J2 ^ [exp (^._i)l exp (-m^._ 4 ) (10) 

j=-M 

137 where 9j = jn/M and the conformally-mapped image function is 

us The argument of f(z) in (111[) is the inverse mapping of z h-s> p, it shows p does not functionally depend on 

139 t, but Weeks' rules-of-thumb for b and k depend on f max . 

«o There are other related methods which use different orthonormal polynomials to represent f(p) inside 

141 the unit circle. Chebyshev polynomials (known as Piessen's method |Piessens(1972)| ) can be used to expand 

142 the f(z) on the real interval [—1, 1]. The Weeks method is moderately easy to program, requiring the use 

143 of Clenshaw recurrence formula to accurately implement Laguerre polynomials. Piessen's method is similar 

144 to implement, with a similar recurrence formula for Chebyshev polynomials. 



145 3.4 Talbot Method 

H6 We can deform the Bromwich contour into a parabola around the negative real axis if f(p) is analytic 
147 in the region between the Bromwich and the deformed Talbot contours |Talbot(1979) . Numerically, j{p) 
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must not overflow as p — > — oo (e.g., in the BEM implementation, the Green's function is the second-kind 
modified Bessel function, which grows exponentially as p — > — oo). Oscillatory /t(p) often have pairs of poles 
near the imaginary p axis; these poles must remain to the left of the deformed contour. 

The Talbot method makes the Bromwich contour integral converge rapidly, since p becomes large and 
negative, making the e pt term in (|4j) very small. A one-parameter "fixed" Talbot method was implemented 
[Abate and Va lk6(2004) ; the Bromwich contour is parametrized as p(0) = r9(cot(9) + i), where < 9 < n, 
and as a rule of thumb r = 2M/(5i ma x). The fixed Talbot method is 



/M0 = ~ 



- N-l 

m e rt + J- *t{eW>)flp{O k ))[l+iC(O k ))} 



fe=l 



(12) 



where ((6 k ) = 6 k + [6 k cot(9 k ) - 1] cot((9 fc ) and 9 k = kn/N [Abate and Valk6(2004)| . Although f(p) doesn't 
depend on t, the free parameter r depends on t m ax- 

Step change ft(p) for non-zero time become very large as p — > — oo, since C [H(t — r)] = e~ Tp /p, where 
Hit — t) is the Heaviside step function centered on time r. This can lead to precision loss, and stability 
or convergence issues with the underlying numerical model, although Mathematica's arbitrary precision 
capabilities have been used to get around this problem [Abate and Valk6(2004)] . 

The fixed Talbot method is very simple to program; [Abate and Valk6(2004) provide a ten-line Math- 
ematica implementation. 



163 3.5 Fourier Series Method 

164 We can manipulate Q into a Fourier transform; first it is expanded into real and imaginary parts (p = 

165 7 + iio), 

7 t 1.00 

/(*) = 2rij t cos M) + «in(orf)] {5R [/(p)] + iQ [/(p)] } idu. 

166 Multiplying out the terms, keeping only the real part due to f(t) symmetry, and halving the integration 

167 range due to symmetry again, leaves 

f(t) = — 3? [/(p)] cosM) - 9 [/(p)] sin(wi) dw. (13) 
77 Jo 

168 When f(t) is real, (|13|) can be represented using the complex form or just its real or imaginary parts. 

169 Although these three representations are equivalent, when evaluating (|13p with the trapezoid rule, the full 



complex form gives the smallest discretization error |Davies(2005)| . The trapezoid rule approximation to 
(|13[) is essentially a discrete Fourier transform, 



/(t,2iV+l) = ^^'K 

k=0 



2N , , 

j I ink \ j nrt 

f I 70 + — I exp I — 



(14) 



172 where 70 = a— log(e)/T, e is the desired relative accuracy (typically 10 -8 in double precision), T is a scaling 

173 parameter (often 2£ max ), and the prime indicates the k = summation term is halved. The p in (|14[) do 

174 not depend on t, but the free parameter T depends on t m ax- 

175 The non-accelerated Fourier series inverse algorithm [TI] is almost useless because it requires thousands 

176 of f(p) evaluations Antia(2002) ( §9.8]. Practical approaches accelerate the convergence of the sum in 1141 

177 Although this is sometimes called a fast-Fourier transform (FFT) method (e.g., Cohen(2007), Chap 4.4]), 
i7s rarely do the number of f(p) evaluations in an accelerated approach justify an FFT approach. The method 
179 implemented uses non-linear double acceleration with Pade approximation and an analytic expression for 

de Hoog et al(1982)| . Although there are several other ways to accelerate the 



the remainder in the series 



Fourier series approach Cohen(2007) , this method is popular and straightforward. Non-linear acceleration 
techniques drastically reduce the required number of function evaluations, but can lead to numerical dis- 
Kano et al(2005),M orales-Casique and Neuman(2009)] . For diffusion, dispersion associated with 



183 persion 



non-linear acceleration is not noticeable. Schapery's, Talbot's, and Weeks' methods are not accelerated in 
a non-linear manner, and therefore may lead to less numerical dispersion, which may be more important 
in wave systems with sharp fronts. 

The creation of the Pade approximation de Hoog et al(1982)| is relatively straightforward in program- 
ming languages that facilitate matrix manipulations (e.g., modern Fortran, Matlab, or NumPy |Oliphant(2007) ). 
There is no dependence on matrix decomposition routines. 



190 3.6 Algorithm Properties Summary 

191 Table [1] summarizes aspects of the five inverse methods. The third column indicates whether p is explicitly 

192 a function of t, the fourth column indicates if the rules-of-thumb used for the optimum parameters depend 

193 on tmax, and the fifth column indicates whether the transform requires complex p and f(p). 

194 For all methods considered here, computational effort to compute f(t) from the vector of f(p) values 

195 was insignificant compared to the effort required to compute the BEM solution used to fill the f(p) vector. 
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Table 1 Algorithmic Summary 



Method 


Limitations on f(p) and fit) 


p(t)? 


P(imax)? 


V 


Stehfest 


no oscillations, no discontinuities in f(t) 


yes 


no 


real 


Schapery 


smoothly varying /(t), / a exists 


no 


no 


real 


Weeks 


none 


no 


yes 


complex 


fixed Talbot 


no high-frequency f(t), f(p — > — oo) exists 


no 


yes 


complex 


Fourier series 


none 


no 


yes 


complex 



This suggests a more complicated method, which allows re-use of f(p) across more values off and converges 
in less evaluations of f(p), would be efficient for Laplace-domain numerical methods. If existing libraries 
or simulations only support real arguments, then the Stehfest, Schapery, or Piessen's methods must be 
used. Complex p methods will pay a slight penalty in computational overhead compared to real-only p 



200 routines. Computing with arbitrary or higher-than-double precision (e.g., Abate and Valk6(2004) ) will 

201 incur a much larger penalty than the change from real to complex double precision. Generally, complex 

202 p methods have better convergence properties than real-only methods. Expansion of f(p) along the real p 

203 axis is separation of non-orthogonal exponentials, while expansion along the imaginary p axis is separation 

204 of oscillatory functions |Lanczos(1988)[ §29]. 



205 4 Numerical Comparison 

206 Four test problems were solved using the BEM for values of p required by each algorithm's rules of thumb. 

207 The test problem domain is a 3 x 2 rectangle, with homogeneous initial conditions and specified potential 

208 at two ends <j>(x = 0) = — 2/t(p), and <j>{x = 3) = 2/t(p), and zero normal flux along the other sides 

209 d(f>/dy(y = {0, 2}) = 0. All plots show the solution computed at a point closer to the x = boundary 

210 (x = 1/3), midway between the insulated boundaries (y = 1). 

211 The first problem computes f(p) using the optimum p at each t (like most inverse Laplace transform 

212 surveys), according to the rules-of-thumb for each method. While this is most accurate, it is very inefficient 

213 - especially when many values of t are required. In the following sections, all methods except Stehfest 

214 use the same f(p) to invert all t. A method's sensitivity to non-optimal free parameters is important in 

215 practical use for Laplace-space numerical approaches. By inverting more than one time with the same set 

216 of Laplace-space function evaluations, large gains in efficiency can be made. The t range used in the plots 
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217 spans three orders of magnitude; it was chosen to show the evolution of potential and substance flux from 

218 initial conditions to steady state. 



219 4.1 Steady Boundary Conditions, Optimum p 

220 The first problem has steady-state boundary conditions. The transient behavior is solely due to evolution 

221 from the zero initial condition, ft(p) = C[H{t)\ = l/p; f(p) has a pole at the origin. All methods performed 

222 equally well with this simple test problem, although the Fourier series method deviates from the finite 

223 difference solutions at larger time. Figure [1] shows the inverted potential and flux using as few evaluations 

224 of f(p) possible, without major deviations from the finite difference benchmark solution. Some trial and 

225 error was needed to use the Schapery method (i.e., further optimization may be possible). 

226 As shown in Figure[2l all the methods performed very well for nine f(p) terms per t but at least 135 f(p) 

227 evaluations are needed total for each method. Schapery's method does the worst in this case, but this may 

228 be improved with further optimization of pj values. The finite-difference approach took at least an order 

229 of magnitude less computational effort for the given accuracy. Making Laplace-space numerical methods 

230 useful alternatives to traditional time-marching approaches, requires improvements to this inefficiency. 

231 4.2 Steady Boundary Conditions, Same p 

232 All methods had more difficulty obtaining accurate results for a wide t range using only one vector of f(p) 

233 (no Stehfest method, since p explicitly depends on t) . Only the last log-cycle of times is inverted accurately 

234 when using nine f(p) (Figure [3]). All the methods - except possibly Schapery's - have a more difficult 
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Fig. 2 Plots of potential and flux through time with five methods for /t(p) = 1/p, using optimum p at each t. 15 X 9 = 135 
total j(p) evaluations are used by each method. Fourier series, Talbot, Stehfest, and Weeks curves are nearly coincident. 




time time 



Fig. 3 Plots of potential and flux through time with four methods for /t(p) = 1/p, same p used across all t. Nine total 
f(p) evaluations are used by each method. 



235 time with the flux at early time (especially the fixed Talbot method). The apparent success of Schapery's 

236 method can be attributed to the expansion of the deviation from steady-state, which in this case decays 

237 exponentially with time. 

238 Figure [3] shows that when increasing to 51 f(p) terms, most convergence problems disappear, except 

239 at small times. Grouping t values by log-cycles and inverting them together using the same /(p) is more 

240 economical than using the optimal p for each t and is still relatively accurate. The results shown in Figure|3] 

241 are nearly as accurate as those shown in Figure [2l but required 1/3 the f(p) model evaluations. 
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242 4.3 Sinusoidal Boundary Conditions, Same p 

243 This problem uses temporally sinusoidal boundary conditions, ft(p) = £(cos4t) = p ^ 16 ■ This boundary 

244 condition violates some assumptions of the inverse transform algorithms (i.e., no steady-state solution, 

245 oscillatory in time), but the behavior is still relatively simple and smooth, with singularities at p = ±4i. 

246 Figure [5] shows the Schapery method fails since there is no / s , but the other methods do well for 19 

247 terms across one t log cycle. Figure [6] shows all methods besides Schapery do well for 51 terms, across 

248 at least two t log cycles. A modified version of (J8]) substituting -Fl 16 for fs/Pj could extend Schapery's 

249 approach to this case, but this solution was not considered here because of its problem specificity. 
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250 4.4 Step-Change Boundary Condition for r > 0, same p 

251 Finally, the same domain was simulated but with step-change boundary conditions at r = 0.08, or ft{p) = 

252 C(H(t — 0.08)) = e _0 08p /p, with singularities at the origin and p = — oo. This function, and those derived 

253 from it (e.g., a pulse or a square wave) are difficult functions to invert accurately, because f(t) is discontin- 

254 uous. Figures [7] and [8] show the Talbot method does not work for t < r in double precision, since ft(p) grows 

255 exponentially as p —> — oo. The Weeks and Schapery methods do worse than the Fourier series approach 

256 (even with N = 51), but their parameters can be optimized further to improve these methods. 

257 Although this step boundary condition could be implemented more accurately by shifting the results 
25s from the first example by t = 0.08, other step-derived time behaviors involving a pulse or square wave 
259 cannot be simplified in this way. 
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Fig. 8 Plots of potential and flux through time with four methods for /t(p) = exp(— 0.08p)/p, same p used across all t. 51 
total /(p) evaluations are used by each method. Weeks' solution is undefined for t < 0.08. 



Table 2 Numerical summary 



Method 


Number of Terms 


Free Parameters 


Implementation 


Stehfest 
Schapery 
Weeks 
fixed Talbot 
Fourier series 


N < decimal precision 
depends on choice of pj 
p — ¥ too slowly as N grows 
p — > — oo quickly as N grows 
p ioo slowly as N grows 


none 

Pj via trial & error 
k & b (very sensitive to b) 
r = (automatic) 
T = 2t max (automatic) 


easiest 
moderate 
moderate 
easy 
most difficult 



260 4.5 Numerical Results Summary 

261 Table [2] summarizes results from numerical testing with these four simple boundary condition time behav- 

262 iors. The second column indicates what limit there is on the number of terms in the approximation and 

263 therefore the accuracy of the method. The size of p required by the Weeks and Fourier series methods grow 

264 much slower than those required by the fixed Talbot method. The third column indicates what parame- 

265 ters are needed to be tuned by the implementer to increase convergence of the method, and whether a 

266 good choice is critical to the success of the method — an automatic method should not require searching 

267 or optimizing parameters to obtain a robust solution. We define robustness as the ability of a solution to 
26a remain useful, even when not at optimality. We prioritize a solution that is good enough and stable over 

269 one that is excellent but catastrophically sensitive to parameter choice. The fourth column indicates the 

270 ease of implementation in modern Fortran, Matlab, or NumPy. The methods could also be implemented in 

271 a variable-precision environment like Mathematica or mpmath [Johansson(2011)] , but this would further 

272 require the BEM model be implemented in such an environment. 
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273 The modest success of the Schapery method is a bit surprising, given its simplicity and use of real p. 

274 The results of the previous section were the product of many iterations of trial and error, this effort was 

275 not included in the implementation effort. A better rule or parametrization of pj might make this method 

276 more widely useful. 

277 The sensitivity of Weeks' method to the parameter choices was also surprising; similarly, the method 



278 could have been improved after some optimization Weideman(1999) , but Weeks' rule of thumb was used for 

279 the parameters. One of the noted advantages of Weeks' method is the need to only compute optimal p once, 
2so then any time can accurately be inverted |Kano et"ld (2005) , Weideman(1999) , Duffy ( 1993) j . When using the 

281 simple rules-of-thumb for the the free parameters, this was not found to be the case. The generalized form 

282 of Weeks' method can include information about behavior of f(p) — > oo (related to behavior as t — > 0) , but 

283 this requires problem-specific knowledge. 

284 The Fourier series method is more robust with respect to non-optimal p values, even though |Duffy (1993)] 

285 cites this as a reason to use Weeks' method over the Fourier series approach. 



286 5 Conclusions 

287 Laplace-space numerical approaches to solve the diffusion equation have several viable alternative inverse 

288 Laplace transform algorithms to choose from. Historically, most Laplace-space solutions to the diffusion 

289 equation have used real-only methods (i.e., Gaver-Stehfest or Schapery). More robust methods require 

290 complex arithmetic and f(p) evaluations, but have the benefits of: 

291 1. handling a broader class of time behaviors (Fourier series method); 

292 2. still being relatively simple to implement (fixed Talbot method) ; 

293 3. only utilizing double-precision complex data types, which are handled natively by modern Fortran, 

294 Matlab, or NumPy, and by common extensions in C++ (Fourier series and Weeks' methods). 

295 Several practical recommendations are made regarding Laplace-space numerical modeling: 

296 1. If many observations are needed across several time log cycles, large gains in efficiency can come from 

297 inverting groups of times with a single f(p) vector (e.g., grouped by log cycle). This complicates the 

298 implementation, but leads to much faster simulations. 
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299 2. If calculating f(p) is very expensive, and some numerical dispersion is allowable (not solving a wave 

300 problem with sharp fronts), then the Fourier series method approach is most economical, and is auto- 

301 matic and robust regarding free-parameter selection. 

302 3. If only a single ft(p) is needed, then it may be worthwhile to optimize free parameters needed by Weeks' 

303 or Piessen's methods, or incorporate information about asymptotic f(p) behavior. Selection of optimum 

304 b is far from automatic, and the Weeks method is not robust for non-optimal free parameters values. 

305 4. If implementation time is a large factor, the fixed Talbot is quite simple to code and was automatic (no 

306 need to select optimum parameters). The fixed Talbot may not work for non-zero step-time behavior 
30? without extended precision. 

308 5. If complex- valued function evaluations are not feasible (e.g., only real matrix or special function libraries 

309 are available) , the Schapery or Piessen's methods are capable of using the same p values to invert different 

310 times, which the Gaver-Stehfest method cannot. 

3u 6. When appropriate, the strategy used by Schapery to expand the deviation from a reference state could 
312 be incorporated as a strategy to improve other algorithms. 
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